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Abstract. In this work we investigate the gravitationally lensed system B1422+231. High-quality VLBI image 
positions, fluxes and shapes as well as an optical HST lens galaxy position are used. First, two simple and smooth 
models for the lens galaxy are applied to fit observed image positions and fluxes; no even remotely acceptable 
model was found. Such models also do not accurately reproduce the image shapes. In order to fit the data 
successfully, mass substructure has to be added to the lens, and its level is estimated. To explore expectations 
about the level of substructure in galaxies and its influence on strong lensing, N-body simulation results of a 
model galaxy are employed. By using the mass distribution of this model galaxy as a lens, synthetic data sets 
of different four image system configurations are generated and simple lens models are again applied to 6t them. 
The difficulties in fitting these lens systems turn out to be similar to the case of some real gravitationally lensed 
systems, thus possibly providing evidence for the presence and strong influence of substructure in the primary 
lens galaxy. 
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1. Introduction 


Gravitational lens systems with multiply imaged quasars 
are an excellent tool for studying the properties of dis¬ 
tant galaxies. In particular, they provide the most ac¬ 
curate mass measures for the lensing galaxy. Besides the 
mass profiles, one can also gain information about evolu¬ 


tion (Kochanek et al. 200C) and extinction laws (Fassnacht 
et al. 1999). Strong lensing is also a very promising and ro¬ 


bust tcol to measure the Hubble constant (Rcfsdai 1964 JT 

The success of this method, however, depends strongly on 
how well the mass model is constrained. 

It turns out that image positions can be f it quite accu- 
rately with simple, smooth elliptical models (Keeton et al. 


1997). Since the number of observational constraints from 

one wants to include the flux 
however, should not be 


image positions is small, 
information. The optical fluxes, 
used, as they mig ht be affected by micro lensing and/or 
dust obscuration (Chang & Refsdal 1979). Radio fluxes, 
on the other hand, can provide further constraints in lens 
modelling. 
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Fitting the fluxes very accurately turns out to be dif¬ 
ficult in many lens systems. Models for MG0414+0534 
(Falco et al. 1997; Keeton et al. 1997), PG1115+080 
(Keeton & Kochanek 1997) and B1422+231 (e.g. 
Kormann et al. 1994Q all show the same failure, namely the 
observed flux ratios are very different from what one would 
expect for the image configurations from a smooth model. 
In particular, the gravitational lens system B1422 +231 
was explored in detail, and as was first mentioned by Mac 


&; Schneider (1998| ), mass substructure in the lens galaxy 


might provide an explanation for the failures in flux mod¬ 
elling. 


A question arises as to whether there is something spe¬ 
cial about these quadruple systems, or does the discrep¬ 
ancy simply arise due to the fact that our smooth mod¬ 
els are oversimplified? In other words, we are asking how 
“well” (for the purpose of strong lensing) can the smooth 
models used to fit the data represent a real galaxy? N- 
body simulation data can provide a realistic description 
of a galaxy mass distribution, thus giving a possibility to 
probe its effect on strong lensing. 


In this paper, we will study the influence of the sub¬ 
structure on the lens system B1422 +231, using VLBI ra - 
dio measurements of the system by Patnaik et al. (1999). 
































2 


M. Bradac et al.: B1422+231:The influence of mass substructure on strong lensing 


In Sec. H we give a description of the lens system and data 
used. The method is outlined in Sect. || and the results on 
fitting the system with smooth mass model are presented. 
In Sect. the model accounting for the substructure is 
presented and in Sect. 5.2 the deconvolved image shape in¬ 
formation is added to the fit. Sect. gives the description 
of the method used to investigate lensing by an N-body 


The lensing galaxy has been observed in the optical; its 
redshift has been determined to be 0.338 and its position 


relative to image B has been measured (Impey et al. 1996). 
The main lens galaxy is a member of a compact group 
with a median projected radius of 35 h kpc and velocity 
dispersion of ~ 550 kms -1 (Kundic et al. 199~^). 


Several groups have tried to model B1422+231 (Hogg 


Simula 

;ed galaxy. We describe how we obtained synthetic 

& Blandford 1994; 

Kormann et al. 1994; 

Keeton et al. 

data 0 : 

four image systems, and the results of fitting such 

1997; 

Mao & Schneider 19981) and all of them have ex- 


systems with lens models are presented. Finally we draw 
some conclusions in Sect. |]. 


This work is an abbreviated version of Bradac (2001 ). 
In the course of writing this paper, several related pa- 
pers on substructure of lens galaxies have been submitted 
(Metcalf fc Madau 2001; Chiba 2001; Dalai fe Kochanek 


2001 ;_ Metcalf & Zhao 2001 Keeton 2001). In Metcalf fe 
Madau (2001 ) and Chiba (2001 ), the authors use semi- 
analytical descriptions to account for mass clumps typical 
for N-body simula tions in a statistical fashion. In addi¬ 
tion, Chiba (2001) tests his prediction on two known sys- 


tems with four images. Dalai fe Kochanek (2001), Metcalf 


& Zha a (2001), and Keeton (2001) predict the proper- 
ties of substructure needed to constrain lens systems and 
compare their results with the CDM predictions. In ad¬ 


perienced difficulties in fitting it. As we used data with 
even more precise image positions one might expect that it 
would become even harder to model the system. However, 
as already pointed out by some authors, the difficulties do 
not lie in fitting the image positions but rather in the flux 
ratios. 

It turns out that simple, smooth models fail to repro¬ 
duce the radio as well as the optical flux ratios of the 
system. While the mismatch in optical data might still 
be due to the microlensing and/or dust obscurations, this 
can probably not explain why such models are not suc¬ 
cessful when fitting radio flux ratios. Mao & Schneider 
(1998|) have proposed a lens model that accounts for the 


dition, Keeton (2001) accounts for the difference between 
radio and optical fluxes by investigating radio and opti¬ 
cal sources of different sizes. The present work is differs 
from the afore mentioned in that slightly different models 
for the lens in B1422+231 system are used and that we 
also include image shape constraints in the fit. Further, 
we investigate four image systems generated by using a 
CDM N-body simulated galaxy, rather than an analytic 
approximation for the statistics of mass-substructure. 

2. The-«*ystei £ y-«f-B4423+251 - 


substructure in the lensing galaxy. They concluded that 
one needs a surface mass density perturbation of the or¬ 
der of 1 % of the critical surface mass density in order to 
change the flux ratios to the observed values. 

In order to succeed in fitting the image flux ratios, one 
needs to consider more sophisticated models. However, 
such models also require the use of additional parame¬ 
ters. Therefore it is very difficult to ensure a constrained 
model that accounts for the substructure using as con¬ 
straints only image positions, flux ratios, and the galaxy 
position. For this reason we also included the axis ratios 
and the orientation angles of the deconvolved i mages as 
additional constraints. They were obtained from Patnaik 


et al. (1999 ) and are given in Table (2| Listed are the ab- 


The gravitational lens system B1422+231 was discovered 
in the course of the J VAS survey (Jordell B ank - VLA 
Astrometic Survey) by Patnaik et al. (1992 ). It consists 
of four image components. The three brightest images A, 
B, and C (as designated by Patnaik et al. 1992 ) are fairly 
collinear. The radio flux ratio between images A and B is 
approximately 0.9, while image C is fainter (flux ratio C 
to B is approximately 0.5). Image D is further away and is 
much fainter than the other images (with flux ratio D:B of 
0.03). The most recent available radio data for the image 
positions and fluxes were obtained from the polarisation 
observations made at 8.4 GHz using the VLBA and the 
100m telescope at Effelsberg from Patnaik et al. (1999) 
and are listed in Table [I]. For each of the components, the 


solute values of ellipticities and the orientation angles of 
the fitted elliptical Gaussians together with uncertainties. 

From the configuration and fluxes of the four images we 
can gain some qualitative constraints on the lens model. 
Image D is much fainter than the rest which can mean 
that it is either highly demagnified or that the other three 
images are highly magnified. Because the position of the 
primary lens galaxy is known, the possibility of image D 
being highly demagnified is ruled out. Namely, image D 
does not lie at a position of high surface mass density 
and thus high demagnification. Images A, B, and C are 
therefore highly magnified. In order to get three highly 
magnified images with a fairly collinear configuration, the 


source has to lie close to (and inside) a cusp (e.g. Schneider 


author j measured positions (relative to the image B) and 
fluxes as well as the deconvolved image shapes. Here and 
through the paper we are using a notation where ( 81 , 62 ) 
are the angular coordinates in the lens plane and (/3i, /? 2 ) 
in the source plane. 8 \ and (3i increase in the negative RA 
direction. 

The radio source of this lens syste m is associated with a 
15.5 mag quasar at a redshift of 3.62 (Patnaik et al. 199^). 


et al. 1992 and references therein). 


For a source position sufficiently close to and inside a 
cusp there exists a relation which states that the sum of 
the fluxes of the outer two magnified images (in our case 
images A and C) is the same as the flux of the m iddle 
image - B (e.g. Mao 1992 ; Schneider fe Weiss 1992 ). This 
rule is strongly violated in the lens system B1422+231. 
Actually, one can get a deviation from this relation if the 
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source is away from the cusp or if there is substructure in 
the system (i.e. fluctuations on a scale smaller than the 
separations between A, B, and C). Also the directions of 
elongation of the images indicate that the source is indeed 
close to a cusp and thus argue in favour of substructure 
in the system (Mao & Schneider 1998). 


Image 

A(9i 
in mas 

A6 2 
in mas 

0RA,Dec 

in mas 

p<i, obs 

mJy 

mJy 

A 

-389.25 

319.98 

0.05 

152 

2 

B 

0.0 

0.0 


164 

2 

C 

333.88 

-747.71 

0.05 

81 

1 

D 

-950.65 

-802.15 

0.05 

5 

0.5 

G 

-717 

-640 

8 




erties. We perform the minimisation in the image plane. 
The x 2 for the positions of the images reads 


Xpos ^ ' 
£= 1 , 3,4 


AG 


, obs 


_ y\^2, mod 


(i) 


2, pOS 


3. Lens modelling 

First we considered two standard gravitational lens models 
since their application to the Patnaik et al. (1999| ) data 
has not yet been discussed in the literature. 

The standard approach to model a strong lens system 
is to define the goodness-of-fit function y 2 , a measure of 
the deviation of the predicted and observed image prop- 


Table 1. Image positio ns with respect to im age B and 
radio fluxes taken from Patnaik et al. (1999| ), where the 
uncertainties of the image positions were estimated to be 
l/20th of the image size in the corresponding direction 
(note that these directions do not coincide with 8 \ and 82 
directions). For simplicity we set the errors to be the same 
in both directions and assign the value of 0.05 mas for all 
images. F l,ohs is the total radio flux density. The position 
of the galaxy (designated by G, measured in the optical) 
was taken from [mpey et al. (1996| ). The uncertainty in 
measuring the galaxy position is much larger than any 
potential misalignment between the optical and radio ref¬ 
erence frames. 


where AG 1 ' obs and AG 1 ' mod are the observed and the mod¬ 
elled position of the *-th image relative to a chosen origin 
(in our case image B, denoted by i = 2), respectively. Note 
that the sum extends only over the images A, C, and D. 

To be more precise, AG 1 ’ mod is obtained as 


AG' 


, mod 


= G l 


mod 


G 


2, mod 


where the vectors G 1 ’ mod and G 2 ’ mod are measured with 
respect to an independent reference point. The same is 
true for AG 1 ' obs , 


2, obs 


obs _ qi, obs _ Q 


which is the observed image position as listed in Table || 
Therefore when the image positions of A, C, and D w.r.t. 
the image B are measured, they all contain the uncer¬ 
tainty of the origin (image B). Thus, the uncertainties of 
the measurements of each of the image positions are cor¬ 
related and we are forced to calculate the x' 2 0 s as i n (jl') • 
We therefore denote 0 + pos to be the uncertainty of the 
measured relative image position. 

Similarly, we define for the galaxy position 


Xgalaxy 


IA< s s - A0™ d | 


o. 


( 2 ) 


gal, pos 


with a ga i t pos being the uncertainty of the relative lens po¬ 
sition. The contribution to the \ 2 -function from the fluxes 
is given by 


Xflux = J2 
2=1 


^^7*2, obs _ jpi, modj' 


( 3 ) 




Table 2. The absolute ellipticity e;| - defined in Eq. (|]) 
and the position angle (measured w.r.t. the $i-axis) of 
the deconvolved image i of the fitted elliptical Gaussians 
calculated from Patnaik et al. (1999| ) data. The uncer¬ 
tainty of the absolute ellipticity o \ t 1 j is determined by 


where F 1 ' mod denotes the flux of the i -th image obtained 
from the model and F 1 ' obs is the observed flux. 

Further we can use deconvolved image shapes as con¬ 
straints of the model. We work in terms of the complex 
ellipticity, 


£i := e, : e 


2i ipi 


taking the uncertainty on the major and minor axis to be 
a tenth of the beam size, which corresponds to 0.1 mas 
and we considered them to be uncorrelated 


Table 3. The summary of parameters used for the lens 
galaxy mass models. 


Model 


Par. 


Description 


Image 

M 


^|e| ,i 

O\p,i 

SIE+SH 

^lens 

OV 

lens position 

Line-of-sight velocity dispersion 

A 

0.70 

143° 

0.07 

5° 


e 

Absolute ellipticity 

B 

0.80 

133° 

0.07 

5° 


<t> 

Position angle of ellipticity 

C 

0.55 

106° 

0.09 

5° 


7i xt ,7l xt 

External shear components 

D 

0.20 

33° 

0.10 

20° 

NIE+SH 

9 C 

Add. to SIE+SH, core radius 
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Each image is thus described by the absolute value of ellip- 
ticity | Ci | and the corresponding position angle ipi (mea¬ 


sured v.r.t. di-axis). From the given data the absolute 
value of ellipticity is calculated as 


uncertainties from Patnaik et al. (1999| ), listed in T able [i~[ 
The optical position of the galaxy was taken from Impey 
et al. (1996|). Although the image positions are very ac- 


ed = 


di - bi 

~bi ’ 




(4) 


where di and bi are major and minor semi axes of the 
fitted elliptical Gaussian. Two additional terms, 


= E 


gobs I _ , e 


mod I 


hi, 


and 


(vf s - <p? od Y 


(5) 


( 6 ) 




have to be added to the ^-function. For simplicity we 
assumed the errors on ellipticity and position angle to be 
uncorrelated. The y 2 function we want to minimise is sim¬ 
ply given by the sum 


2 2,2 i 2 i 2 i 2 

Xtot Xpos Xgalaxy "f - Xflux d - Xe d“ Xt^ 


(7) 


In order to obtain the values of A 9 h mod from the 


lens equation (see e.g. (Schneider et al. 1992|), we used 
the Numerical Recipes MNEWT routine from Press et al. 
(1992[)]for solving a set of two non-linear equations. The 
y 2 -function was then minimised with respect to the model 
parameters using POWE LL, a multi-dimen sional minimisa¬ 
tion routine, also from Press et al. (1992|) . 

One of the consideration one faces when doing lens 
modelling is the uniqueness of the solution. Minimisation 
methods such as POWELL do not give an answer as to 
whether a minimum found is a local or global one. One 
can partly solve this by using simulated annealing rou- 


curate (of the order of 50 /zarcsec), we have no difficulties 
fitting them, and so the y 2 contribution from the image 
positions drops to zero (see Table ||). However, as already 
pointed out in previous works on B1422+231, the model 
completely fails in predicting the image fluxes. In particu¬ 
lar image A is predicted too dim (the modelled flux ratio 
A:B turned out to be 0.80, much below the measured value 
of 0.93). We have also tried to model the system with a 
NIE+SH model; however, the y 2 did not improve signifi¬ 
cantly. 

Next we disregarded the flux of image A as a con¬ 
straint, trying to see whether we can get a good fit for the 
rest of the images. Reducing the number of constraints by 
one, we are left with one degree of freedom. The results of 
the fitting are given in Table [|. The fit is still not satis¬ 
factory, thus indicating that changes of magnifications are 
needed not only for image A. 

The choice of these two particular macromodels was 
due to their simplicity and analytic solutions for the de¬ 
flection angle. Of course these are not the only models one 
can use, and one might hope to explain the discrepant flux 
ratios by using another family of models. In order to effec¬ 
tively change the flux ratios, one has to change the magni¬ 
fication gradient. The simplest way of changing a gradient 
is to change the power-law behaviour of the potential. 

Although the success of such modelling is very unlikely 

we have tried to 


(as shown by Mao & Schneider IE 
use the SIE+SH model, adding an additional black hole 
with an Einstein radius fixed to 10 -2 of the main galaxy 
Einstein radius, centred on the lens galaxy position. The 
corresponding mass of the black hole is ~ 1O 8 M 0 , in agree¬ 
ment with the correlation of the black hole mass with the 


tines, such as AMOEBSA from Press et al. (1992|). However, 
these routines are much more CPU-time consuming, and 
they are also not completely foolproof in finding an abso¬ 
lute minimum. 

As an additional check, to minimise the chance of ob¬ 
taining a secondary minimum as a result, we have run 
POWELL for several different initial conditions. Unless the 
global minimum is extremely narrow in parameter space 
(which is rather unlikely) we can be confident of obtaining 
the global minimum. 


velocity dispersion of its host bulge (Ferrarese fc Merritt 


2000| ). Such a model is plausible, however we found no sig¬ 
nificant improvement of the fit (see Table f| for results). 
Also, considering the black hole mass as a parameter did 
not significantly improve the outcome. 

Although other smooth models can still be investi¬ 
gated, it seems unlikely that another smooth model can 
explain all four image fluxes simultaneously. Even when 
one disregards the flux of image A, a smooth macromodel 
seems to be incapable of explaining the remaining flux 
ratios. 


4. Modelling B1422+231 with smooth models 

We used a singular isothermal ellipsoid with external 


shear from Kormann et al. (1994) (hereafter SIE+SH) and 
a non-singular isothermal ellipsoid model with external 
shear (NIE+SH) from Keeton fe Kochanek (199~§| ) to fit 
the image positions and fluxes of B1422+231. The expla¬ 
nation of the model parameters for the models we used 
are given in Table [| 

We have applied the fitting procedure described above 
to the radio data, using image positions, fluxes, and their 


5. Models with substructure 

In the previous section it turned out that A:B flux ra¬ 
tio causes the biggest difficulty in fitting the B1422+231. 
Since the radio and optical flux ratios are very different, 


one is tempted to exclude it from the measure (see Mac 
|fe Schneider 199S , Chiba 2001). 


However, one can also try to deal with this problem in 
another way. Adding a small perturber at the same angu¬ 
lar diameter distance as the primary lens and at approxi¬ 
mately the same position as image A can change the flux 
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ratio A:B substantially. On the other hand, calculations 
show (Mao & Schneider 1998) that such a perturber does 
not affect the positions of any of the images appreciably. 
Furthermore, a small perturber can also change the flux 
ratio of the other images slightly and this might help to 
improve the results from the previous section. 

We model the perturber as a non-singular isothermal 
sphere. This gives two additional parameters (two per¬ 
turber positions, since we keep the line-of-sight velocity 
dispersion and core radius of the perturber fixed) to the 
macro (SIE+SH) model. The choice of the perturber being 
modelled as non-singular is due to the fact that a singular 
isothermal sphere (with the same Einstein radius) is more 
likely to give rise to additional (observable) images. 

We are aware of the fact that the choice of this partic¬ 
ular model for the substructure is oversimplified in many 
ways. However, we are not trying to constrain the nature 
of substructure in this case; which is impossible due to 
the number of constraints available. In the following we 
investigate whether the peculiarity of the fluxes can be 
explained by a small satellite galaxy in the neighbourhood 
of one (or two) of the images. 


5.1. Modelling B1422+231 using the substructure 
model 

For modelling B1422+231 we fixed the line-of-sight veloc¬ 
ity dispersion of the NIS perturber to onis = 10 kms -1 
(equivalent to an Einstein radius of approximately 2 mas) 
and its core radius to 9 C = 20 mas. A perturber with 
these properties does not affect the image positions signif¬ 
icantly; on the other hand it can substantially change the 
magnification at the position of one of the images. 

When fixing the core radius one has to be aware that 
not only an SIS, but also an NIS lens with small enough 
core radius might give additional images. Therefore, an 
NIS perturber should have a core radius much bigger than 
the Einstein radius, in order not to produce additional im¬ 
ages. We have checked that indeed no additional observ¬ 
able images are predicted by the model. 

For this purpose we define the function / on a grid of 
points 9j 

f = \P(9 j )-f3 s \ 2 , (8) 

where [3(0 j) is the calculated source position correspond¬ 
ing to the point Oj. The position (3 S is defined as the aver¬ 
age position for the source predicted by the four observed 
image positions 

i 4 

A = i E^(^’ obs )- (9) 

The function / vanishes only around the positions where 
images are observed, provided that the grid resolution is 
high enough. In our case the grid resolution was 0.05 mas 
(small compared to the Einstein radius of the perturbing 
galaxy, 0e,nis = 3 mas). We indeed found only four such 


regions, and they correspond to the four observed images. 

The resulting model has 12 parameters, which leaves us 
0 degrees of freedom. The y 2 has decreased by a factor of 
more than 20 compared with the SIE+SH model; however, 
since we have zero degrees of freedom we expect \ 2 t° 
vanish if the model is realistic (and if the x 2 technique 
is an adequate method). The family of models considered 
thus does not seem to be adequate for the description of 
the galaxy in B1422+231 lens system. 

We again see that the flux ratio of images A:B is not 
the only problem when dealing with fluxes. Adding a small 
perturber close to image A allows us to adjust the mod¬ 
elled flux of A such, that it does not give any contribution 
to the % 2 -function. Still the remaining image fluxes are 
not fit perfectly, leading to a possible conclusion that all 
images are affected by the mass substructure. 


5.2. Using image shapes as constraints 


In order to increase the number of degrees of freedom we 
are dealing with, further constraints have to be added. We 
therefore used the deconvolved image shapes (see Table |j) . 
Due to the problems in determining the image shapes we 
describe later, we will use them as constraints only in this 
section. 

The \ 2 minimisation was done according to the pre¬ 
vious section, and the results are presented in Table 
We have decreased the value of the core radius of the per¬ 
turber in order to get higher magnification gradients, since 
large changes in image shapes as predicted by a smooth 
model was needed. Again, no additional observable images 
are produced by the perturber(s). 

In the Patnaik et al. (19~99| ) paper the uncertainties on 
the image shapes are not listed. The image shapes are ob¬ 
tained by fitting Gaussian profiles to the map, and then 
deconvolved using the known beam-shape. The uncertain¬ 
ties are therefore just a rough estimate, since one can not 
quantitatively account for the error of such fitting. In fact, 
all the images exhibit non-Gaussian features a Gaussian 
model is an oversimplification for the image shape descrip¬ 
tion. 


Table 4. Result of modelling B1422+231 radio positions 
and flux data with (i) an SIE model with external shear 
(ii) SIE+SH without the flux of image A as a constraint 
and (iii) SIE+SH model with an additional black hole at 
the galaxy position without the flux of image A as a con¬ 
straint. The black hole Einstein radius in model (iii) was 
fixed to 10~ 2 of the main galaxy Einstein radius. The pa¬ 
rameters of the best fitting model are listed in Table || 


Model 

Adof 

Xtot 

Xpos,rel 

”1” Xgal,rel 

+ 

Xflux 

(i) 

2 

129 

= 0 

+ 

18 

+ 

111 

(ii) 

1 

9.8 

= 0 

+ 

0.5 

+ 

9.3 

(iii) 

1 

7.5 

= 0 

+ 

0.2 

+ 

7.3 
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The resulting \ 2 for the minimisation was 19, having 6 
degrees of freedom (the probability of obtaining a value for 
X 2 bigger than 19 is 0.0042). We further try to fit the data 
with two perturbers in the system. If we put two equal 
perturbers into the system (i.e. both with fixed 9 C = 1 mas 
and (Tnis = 10 krn/s), we get a \ 2 of 15 (see Table[|) with 
4 degrees of freedom. The probability of obtaining a % 2 
larger than that value is now 0.0047. Such a reduction of 
X 2 is apparently not a significant improvement of the fit 
(compared to the model with a single perturber). 

Just for comparison, we also include the deconvolved 
image shapes in the fit with SIE+SH model. The result¬ 
ing x 2 is 140 for 8 degrees of freedom (the probability 
of obtaining a value for \ 2 larger than that is now only 
10~ 28 ). We essentially get the same model as when only 
fitting image positions and fluxes (see Table ||). We see 
that the fluxes provide much stronger constraints (their 
uncertainty is much smaller) than the image shapes. 

Apparently, models with substructure yield signifi¬ 
cantly better fits than the ones without. However, we 
should stress that, strictly speaking, the appropriate sta¬ 
tistical treatment can not be easily performed because the 
model parameters enter the y 2 function in a non-linear 
way. As a result, the % 2 function does not behave as a 
chi-square random variable with the appropriate number 
of degrees of freedom. One can see that already from the 
fact that the \ 2 function did not vanish when using a 
model with zero degrees of freedom. Hence, the probabil¬ 
ities we quoted above have to be taken with care; still 
they can be used to compare the performance of different 
models. 

It is clear that it is difficult to simultaneously stretch 
and rotate the images with one (or two) perturber(s). The 
macro-model is not very successful in predicting both com¬ 
ponents of the image ellipticities and the fluxes, and there¬ 
fore corrections are needed in the case of all four images. 
We can safely assume that the inclusion of further sub¬ 
clumps in the model would eventually lead to a “perfect” 


Table 5. Result of modelling B1422+231 radio positions, 
flux and in the case of (v), (vi) and (vii) image shape 
data. The models we used were an SIE model with exter¬ 
nal shear and an additional perturbing NIS galaxy (iv), 
an SIE model with external shear (v), SIE+SH and one 
additional perturbing NIS galaxy (vi), and SIE+SH with 
two perturbing galaxies (vii). The core radius and <tnis 
of the perturbing galaxy(ies) were fixed to the values 
0 C = 20 mas and ctnis = 10 km/s for model (iv) and 
9 C = 1 mas and ctnis = 10 km/s for the other models (see 
text). The parameters of the best fitting models are listed 
in Table 


Model 

IVdof 

Xtot 

Xpos,rel 

+ 

Xgal,rel 

Xflux 

+ Xi + X% 

(iv) 

0 

5.6 

= 0.0 

+ 

0.8 

+ 4.8 


(v) 

8 

140 

= 0 

+ 

18 

+ 111 

+ 4+7 

(vi) 

6 

19 

= 0 

+ 

0 

+ 6 

+ 9 + 4 

(vii) 

4 

15 

= 0 

+ 

0 

+ 4 

+ 7 + 4 


fit with the observed data. In particular three more sub¬ 
clumps close to the images B, C, and D would yield a 
significant improvement to the y 2 . 


6. Strong lensing by an N-body simulated galaxy 

A question that arises from model fitting of B1422+231 is 
whether such behaviour is seen by the N-body simulated 
galaxy and therefore generic of a typical galaxy lens. We 
used the cosmological N-body sim ulation data including 
gasdynamics and star formation of Steinmetz & Navarro 


(200l|) for this purpose. 


The simulations were performed using GRAPESPH, 
a code that combines the hardware N-body integrator 
GRAP E with the Smo oth Particle Hydrodynamics tech¬ 
nique ( [Steinmetz 1996|) . GRAPESPH is fully Lagrangian 
and optimally suited to study the formation of highly non¬ 
linear systems in a cosmological context. The version used 
here includes the self-gravity of gas, stars, and dark mat¬ 
ter components, a three-dimensional treatment of the hy¬ 
drodynamics of the gas, Compton and radiative cooling 
(assuming primordial abundances), the effects of a photo- 
ionizing UV background, and a simple recipe for trans¬ 
forming gas into stars. 

The original simulated field is located at z = 0.2 and 
contains approximately 300000 particles. The simulation 
is contained within a sphere of diameter 32 Mpc which 
is split into a high resolution sphere of diameter 2.5 Mpc 
centred around a galaxy and an outer low resolution shell. 
Gas dynamics and star formation is restricted to the high 
resolution sphere (280000 particles, 92000 of which dark 
matter), while the 34000 dark matter particles of the low 


Table 6. Resulting parameters from best fitting models: 
(i) SIE+SH, (ii) SIE+SH, (iii) SIE+SH+BH, both with¬ 
out the flux of image A as a constraint, (iv) SIE+SH+NIS, 
(v) SIE+SH, (vi) SIE+SH+NIS, (vii) SIE+SH+2NIS. In 
models (v), (vi) and (vii) we use the shapes of the images 
as additional constraints. #nis is the position of the per- 
turber(s) w.r.t. the image indicated in the superscript. The 
parameters of the best fitting source shape were |e s | =0.14 
and ip s = 60° for the model (iii), for the model (iv) 
|e s | = 0.04, ip s = 30°, and for the model (v) |e s | = 0.10, 
( p s = 20°. The position angles are measured w.r.t. 01- 
axis. The resulting line-of-sight velocity dispersion er v was 
190 kins -1 in all four cases. 



^lens 

Cgal 


ext 

7i 

ext 

72 

0NIS 


(mas, 

mas) 





(mas, mas) 

(i) 

(-744, 

-659) 

0.19 

34° 

-0.04 

-0.16 


(ii) 

(-721, 

-643) 

0.14 

33° 

-0.05 

-0.17 


(hi) 

(-720, 

-643) 

0.15 

33° 

-0.05 

-0.18 


(iv) 

(-722, 

-646) 

0.13 

32° 

-0.05 

-0.18 

(41, 36) A 

(v) 

(-744, 

-659) 

0.19 

34° 

-0.04 

-0.16 


(vi) 

(-718, 

-643) 

0.12 

32° 

-0.05 

-0.18 

(53, 47) a 

(vii) 

(-719, 

-641) 

0.13 

32° 

-0.05 

-0.18 

(—18, 7) a 
(-5,4)° 
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Fig. 1. The cut-out of the surface mass density map of the 
simulated galaxy. The mass distribution resulting from the 
cosmological N-body simulation (see text) was smoothed 
using convolution with a Gaussian kernel characterised by 
a standard deviation a ~ 0.8 kpc ~ 0.2 arcsec. This map 
was then evaluated on 2048 x 2048 grid points (~ 160 x 
160kpc) and the surface mass density was calculated. Test 
particles, used in the N-body simulation to account for 
the large-scale structure, have been removed here. The 
contours correspond to the values of k = 0.8; 1.6; 2.4; 3.2. 
The dark regions represent the regions of high n. The units 
on the axes are arcseconds; one arcsecond in the lens plane 
corresponds to approximately 4 kpc. 


resolution sphere sample the large scale matter distribu¬ 
tion in order to appropriately reproduce the large scale 
tidal fi elds (see Navarro & Steinmetz 1997 and Steinmet^ 


& Navarro 200C for details on this simulation technique). 


The simulation was performed in a ACDM cosmology 
(fl 0 = 0.3, 0\ = 0.7, fl b = 0.019 /h 2 , erg = 0.9). It has a 
mass resolution of 1.26 x 1O 7 M 0 and a spatial resolution 
of 0.5kpc. A realistic resolution scale for an identified sub¬ 
structure is typically assumed to be ~ 40 particles which 
corresponds to 5 x 1O 8 M 0 . The quoted mass resolution 
holds for gas/stars. The high resolution dark matter par¬ 
ticles are about a factor of 6 (= fio/^b) more massive. 

From the original simulated field we took a cut-out 
map of size ~ 160 x 160 kpc that is centred on a sin¬ 
gle galaxy. This area lies well within the high resolu¬ 
tion sphere and is void of any massive intruder particles 
from the low resolution shell. The resulting mass distri¬ 
bution was smoothed using convolution with a Gaussian 
kernel characterised by a standard deviation of a ~ 
0.8 kpc ~ 0.2 arcsec.^ This map was then evaluated on 
2048 x 2048 grid points. The final map contains infor¬ 
mation about approximately 130000 particles with a total 


1 For the distance calculations through the paper we assumed 

an Einstein-de-Sitter Universe and the Hubble constant Ho = 

65 kms -1 Mpc -1 . 


mass of 3.0 x 10 12 M 0 . The surface mass density k was 
calculated for every grid point. We chose the redshift of 
the source to be z = 3. A part of the cut-out map can be 
seen in Fig. [I], 

The lens properties are calculated on a grid of 2048 x 
2048 points. The Poisson equation 


V 1 2 ip(d) =2k(0) . 


( 10 ) 


is solved (on the grid) in Fourier space by the 
FFT (Fast Fourier Transformation) method. 
This method is incorporated in the routine 
KAPPA2STUFF from the IMCAT software of Nick Kaiser 


(http://www.ifa.hawaii.edu/~kaiser) that we used. 
It takes a grid map of n(9j) as an input and returns 
the values of deflection angle and complex shear (in 
cc-space), along with other quantities. From these data 
we calculated the Jacobi matrix for each grid point. 

The simulated galaxy is a field galaxy. Therefore we 
add two external shear components to the Jacobi matrix 
(evaluated at each grid point) in order to make it similar 
to the galaxy in B1422+231. The shear components were 
taken to be the same as the ones obtained from the best 
fitting SIE+SH model 


7® xt = —0.04 , 72 Xt = —0.16 . 


The external shear accounts for the effect of the neigh¬ 
bouring galaxies of the compact group, which are not 
present in the simulation. 



Fig. 2. The magnification map of the simulated galaxy 
calculated using the KAPPA2STUFF routine from Nick 
Kaiser’s software IMCAT. External shear is added in the 
evaluation of the magnification map for to account for 
neighbouring galaxies (see text). Lighter regions represent 
high magnifications. The units on the axes are arcseconds; 
one arcsecond in the lens plane corresponds to approxi¬ 
mately 4 kpc. 
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Fig. 3. The caustic obtained by mapping the points of 
high magnification onto the source plane. We took \fj,\ > 30 
for the inner caustic, for the outer one we additionally 
picked the points of \fi\ > 0.5 (due to nummerical effects - 
see text) from the central part of the magnification map. 
The units on the axes are arcseconds; one arcsecond in the 
source plane corresponds to approximately 6 kpc for the 
source at z = 3. 

Fig. 1 shows the magnification map of the surface mass 
density (given in Fig. [l|) with additional external shear. 
One can clearly see the outer critical curve (white curve), 
while only the traces of the inner critical curve are visible 
(little circle inside the black region). The reason why we 
can see the outer critical curve so much better than the 
inner one is the following. At the centre Q c of the galaxy, 
the surface mass density is very high and the determinant 
of the Jacobi matrix can be approximated by detA ~ 
k( 6 c ) 2 . In our particular case, k (0 c ) « 50, and since at 
the critical curve det A = 0, we see that the determinant 
has to decrease from 2500 to 0 in a region of 0.4 arcsec. 
The transition is therefore very steep and we have a very 
good chance to miss the maximum value of magnification 
there, since the resolution is not high enough. At the outer 
critical curve, the change is slower and we can clearly see 
the points of high magnification. 

In order to generate a similar image configuration as 
the one in B1422+231, one considers the caustic curve. 
This can be done by simply mapping the points of high 
magnification onto the source plane. Such a map is pre¬ 
sented in Fig. |j. 

On first sight the caustic structure of the N-body sim¬ 
ulated galaxy looks the same as e.g. the caustic of the 
smooth NIE model with a small core radius. However, if 
we look only at the inner, asteroid caustic we can see that 
it is not completely “smooth”. With the help of bilinear 
interpolation we recalculated the magnification map on a 


Fig. 4. The caustic curves obtained by first interpolating 
the magnification map on a refined grid using bilinear in¬ 
terpolation (increasing the number of points in the region 
of interest by a factor of 25) and mapping the points of 
high magnification |/i| > 30 (Fig. ||) to the source plane. 
The units on the axes are arcseconds; one arcsecond in 
the source plane corresponds to approximately 6 kpc for 
the source at z = 3. Two marks correspond to the source 
position of data set 1 (cross) and 11 (star) - see Sect. |0.l[ 
The source positions of other data sets are located on a 
line connecting them. 


refined grid (increasing the number of points being evalu¬ 
ated by 5 x 5) and the corresponding caustic for such a 
grid is shown in Fig. |]. The caustic structure is much more 
complicated than in the case of a smooth model; in addi¬ 
tion to folds and cusps we also have swallowtails formed 
(see [Schneider et al. 1992 for more details on catastrophe 
theory). 


6.1. Generating synthetic data 

We select a source position (3 S such as to lie inside the 
asteroid caustic and close to the cusp, trying to chose a 
position for which we would get similar flux ratios as in 
the case of B1422-)-231. In total we considered 11 different 
source positions (see Fig. Q). 

For each of them we first determine approximate im¬ 
age positions using the method described in Sect. 5.1. In 


order to get exact image positions we use the root find¬ 
ing method MNEWT again, for which we need the deflection 
angle to be continuous inside the region where we look 
for images. In our case the deflection angle is defined only 
on a 2048 x 2048 grid. We perform bilinear interpolation 
between grid points. Having the image positions, we per- 
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formed bilinear interpolation of the magnification map in 
order to get more exact magnification factors. 


6.2. Fitting the synthetic data 

For the \ 2 -fitting method according to ( 0 ) we need to de¬ 
termine the uncertainty on the image positions and fluxes. 
Since we use interpolation for the MNEWT method we do not 
have a real estimate for the errors. One can, for example, 
set the errors to the same (relative) values as the uncer¬ 
tainties on the observed radio positions in B1422+231. For 
the typical scales we are using here (i.e. the distance B to 
D is approximately 3.5 arcsec) this would mean an un¬ 
certainty of much less than a distance between two grid 
points. However, such a small error estimate is not realis¬ 
tic; due to the finite grid we estimate the image position 
uncertainties to be the distance between two grid points 
(which is a generous estimate; we use bilinear interpola¬ 
tion so the uncertainty is probably lower). The flux ratio 
errors were then set to be approximately 2 % for images 
A, B and C and 5 % for image D. These uncertainties are 
set to be the same as in the case of observed radio fluxes 
in B1422+231. The galaxy position error is set twice as 
big as that for the image positions. 

The fitting procedure is performed in the same way 
as for the B1422+231 data. Again, image B is taken as 
a reference and the y 2 -function is evaluated according to 
( 0 )- We try to fit the positions and fluxes with SIE+SH 
and SIE+SH+NIS models. The flux ratios of the 11 sets 
of synthetic data, together with the results of the min¬ 
imisation are presented in table [l]. We experience similar 
problems fitting fluxes as before; the y 2 -function value is 
high for all 11 data sets. 

What might be surprising is the fact that we do not re¬ 
cover some properties of the lensing system. We know that 
the primary lens is fairly circular (one can not see that 
from fig. since there, external shear is already added) 
and we know the values of the external shear components. 
What we also know a priori are the magnification factors 
for the images. These values were not recovered with high 
accuracy in model fitting. 

For a smooth model and a source close to the cusp, 
the flux relation described before holds. We see that the 
fluxes violate this rule in all configurations we used. As we 
mentioned before this relation can only be violated when 
the source is away from the cusp or if there is substructure 
in the system. Since here we know the source position, the 
N-body lensing results show that the substructure we have 
in this particular simulated galaxy is indeed responsible 
for the observed deviation. 


6.3. Discussion of the results from N-body lensing 

An important question is whether the N-body simula¬ 
tion galaxy we are using is a good representation of a 
real galaxy for the purpose of lensing. If the resolution is 


not high enough, an N-body simulated galaxy might show 
more substructure than a real galaxy has. 

In order to obtain the surface mass density map rep¬ 
resentative of lensing and to try to make sure that the 
substructure we see is not of numerical origin we used a 
smoothing length for particles of a = 0.8 kpc ~ 0.2 arcsec. 
The individual mass clumps we see in the corresponding 
surface mass density map (Fig. [l]) therefore contain well 
above 100 particles. As we mentioned before, a realistic 
resolution scale for an identified substructure in the sim¬ 
ulation corresponds to ~ 40 particles. 

The regions that are interesting for multiple image for¬ 
mation typically have n values of order unity. In such re¬ 
gions we have approximately 300 particles per smoothing 
area, which, if the noise were Poissonian, would result in 
an error of about 1/y/N ~ 0.06. However, studies of errors 
of N-body simulations using smoothed particle hydrody¬ 
namics ( Monaghan 1992| ; |Niedereiter 1978 ) show that the 
errors are significantly smaller, and behave as oc log N/N. 

Any deviation of k due to substructure larger than the 
estimated error of the N-body simulations tells us that we 
are dealing with “physical” substructure, i.e. not of numer¬ 
ical origin. Our surface mass density maps show deviations 
well above the estimated Poisson noise, and therefore also 
above numerical noise. The results shown in Table [l] are 
affected by the “physical” substructure as well as by the 
noise. Although the expected noise is significantly smaller 
than believable substructure observed in the N-body sim¬ 
ulations, it is hard to precisely estimate the relative con¬ 
tributions in the results shown here. 

There are indications that the level of substructure as 
obtained from simulations can influence lensing phenom¬ 
ena a lot. In particular, the synthetic fluxes we obtained 
deviate strongly from those predicted by smooth mod¬ 
els. This particular example of a simulated galaxy can of 
course not give us the answers to the aforementioned ques- 


Table 7. The flux ratios of image A, C and D w.r.t. image 
B of 11 data sets with different image positions (3 S (see also 
Fig. §. Listed are also the resulting \ 2 values for fitting 
image positions and fluxes with SIE+SH model {x\) and 
SIE+SH+NIS model (Xi)- The core radius and unis of 
the perturbing galaxy were fixed to the values 9 C = 1 mas 
and ctnis = 15 km/s. 


Data set 

Fab 

Fcb 

Fdb 

xi 

X2 

1 

1.04 

0.80 

0.220 

120 

2.0 

2 

1.29 

1.20 

0.193 

960 

120 

3 

1.43 

1.40 

0.167 

2100 

190 

4 

1.20 

0.79 

0.142 

660 

5.1 

5 

1.06 

0.59 

0.116 

350 

3.9 

6 

1.05 

0.42 

0.089 

150 

81 

7 

1.06 

0.31 

0.064 

450 

110 

8 

0.60 

0.29 

0.047 

93 

40 

9 

0.51 

0.46 

0.051 

400 

180 

10 

0.59 

0.52 

0.064 

86 

9.3 

11 

0.62 

0.66 

0.070 

240 

36 
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tions. In order to draw stronger conclusions, one would 
have to investigate many different realisations of N-body 
simulated galaxies and in addition use higher resolution 
simulations (currently unavailable). A statistical analysis 
to investigate the strong lensing properties could then be 
made. This is, however, beyond the scope of this paper 
and will be considered in a future work. 

7. Conclusion 


the fluxes of more than one image are probably affected 
by the lens substructure. With the conclusions one can 
draw from lens modelling we can not give a direct answer 
about the properties of mass substructure in B1422+231. 

One should therefore avoid using the flux constraints 
directly; they should, rather, be treated in statistical man¬ 
ner, e.g. in a way suggested by Mao & Schneider (1998). 
Fortunately, the perturbations on the scales we are deal¬ 
ing with here do not influence the image positions si gnifi¬ 
cantly, and play even less of a role for the time delay ( Mao 


In this work we have investigated the influence of sub- fe Schneider 1998|). Strong lensing thus remains one of the 


structure in the gravitationally lensed system B1422+231. 
While it is intuitively clear that a lens galaxy is not a 
smooth entity, we have tried to investigate how deviation 
from a smooth model can influence lensing phenomena, 
especially the image flux ratios. 

We have used two different smooth models for the lens¬ 
ing galaxy (SIE+SH and NIE+SH), and both failed very 
badly in fitting the image fluxes (we got y 2 = 130 with 2 
degrees of freedom). The use of models with substructure 
requires additional observational constraints. Therefore, 
we used deconvolved image shapes as constraints. We get a 
significant improvement of the fit compared to the smooth 
model. However, the way the substructure is introduced 
is oversimplified, thus we should not be surprised that 
the resulting % 2 is still high. For the model with a sin¬ 
gle perturber we got y 2 = 19 for 6 degrees of freedom, 
and with two perturbers we had y 2 = 15 for 4, while 
the model without substructure (where deconvolved im¬ 
age shapes were included) gives % 2 = 140 for 8 degrees of 
freedom. 

Up to now we have not considered the possibil- 
ity that microlensing plays a role for the radio fluxes. 
Koopmans fe de Bruyn (2000 ) claim that they have de¬ 
tected microlensing in the multiply-imaged radio source 
B1600+434. Microlensing is a very tempting explanation 
for difficulties in fitting the fluxes for it can also explain 
why t he 8.4 GHz A:B flux ratio has change d from 0.97 in 
1991 ( |Patnaik et al. 1992 ) to 0.93 in 1997 ( Patnaik et al. 


1999| ). This has a consequence that again speaks in favour 


of substructure, since the presence of radio microlensing 
indicates that there is a significant number of compact 
objects in the lens galaxy halo. 

N-body simulation data of a model galaxy provides a 
test for the influence of mass-substructure in strong grav¬ 
itational lensing. When we generated data of four image 
systems with the simulated galaxy we again experienced 
difficulties. We have tried to fit image positions and fluxes 
and failed to obtaining a model that fits well. From these 
experiments we can conclude that the level of substructure 
obtained from this particular N-body simulated galaxy 
can cause the same difficulties as experienced in some of 
the real gravitationally lensed systems. 

In order to obtain stronger conclusions one would have 
to investigate more realisations of simulated galaxies, also 
at different redshifts. However, the N-body simulation re¬ 
sults indicate that substructure plays an important role 
in strong lensing. Also, modelling B1422+231 shows that 


best tools to constrain the Hubble constant. 
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